Load data
## Load single-cell experiment
sce <- readRDS(opt$sce)
sce.unstim <- sce[["unstim"]]
sce.stim <- sce[["stim"]]
rm(sce) # free-up the memory
## Load cell type annotation
types <- read.csv(opt$types) %>%
mutate(uniqueBarcode = as.character(uniqueBarcode)) %>%
mutate(celltype_broad = ifelse(celltype %in% c("CD16- NK-cell", "CD16+ NK-cell"), "NK-Cells", celltype)) %>%
mutate(celltype_broad = ifelse(celltype %in% c("T-reg", "gd T-cell", "naive CD4 T-cell", "naive CD8 T-cell", "CD8 T-emra", "CD8 effector memory T-cell",
"CD8 central memory T-cell", "CD7+ memory CD4 T-emra", "CD7+ effector memory CD4 T-cell",
"CD7+ central memory CD4 T-cell", "CD7- memory CD4 T-emra", "CD7- effector memory CD4 T-cell",
"CD7- central memory CD4 T-cell"), "T-Cells", celltype_broad)) %>%
mutate(celltype_broad = ifelse(celltype %in% c("naive B-cell", "memory B-cell"), "B-Cells", celltype_broad)) %>%
mutate(celltype = ifelse(celltype %in% c("CD7- effector memory CD4 T-cell", "CD7+ effector memory CD4 T-cell"), "CD4 EM T-cell", celltype),
celltype = ifelse(celltype %in% c("CD7- central memory CD4 T-cell", "CD7+ central memory CD4 T-cell"), "CD4 CM T-cell", celltype),
celltype = ifelse(celltype %in% c("CD8 central memory T-cell"), "CD8 CM T-cell", celltype),
celltype = ifelse(celltype %in% c("CD8 effector memory T-cell"), "CD8 EM T-cell", celltype),
celltype = ifelse(celltype %in% c("CD7- effector memory CD4 T-cell", "CD7+ effector memory CD4 T-cell"), "CD4 EM T-cell", celltype),
celltype = ifelse(celltype %in% c("CD7- memory CD4 T-emra", "CD7+ memory CD4 T-emra"), "CD4 T-emra", celltype))
## Drug screen data
drugList <- readRDS(opt$drugscreen)
## Define drug combinations
drugComb <- c("Birinapant|Necrostatin-1 (25)", "Birinapant|QVD-Oph (25)", "Birinapant|Necrostatin-1 (25)|QVD-Oph (25)",
"Birinapant|Necrostatin-1 (12.5)", "Birinapant|QVD-Oph (12.5)", "Birinapant|Necrostatin-1 (12.5)|QVD-Oph (12.5)",
"GDC-0152|Necrostatin-1 (25)", "GDC-0152|QVD-Oph (25)", "GDC-0152|Necrostatin-1 (25)|QVD-Oph (25)", "GDC-0152|Necrostatin-1 (12.5)",
"GDC-0152|QVD-Oph (12.5)", "GDC-0152|Necrostatin-1 (12.5)|QVD-Oph (12.5)", "Birinapant|Ipatasertib (2)",
"Birinapant|Ruxolitinib (2)", "Birinapant|Dacinostat (0.04)", "Birinapant|Venetoclax (0.04)",
"Birinapant|Bafilomycin_A1 (0.08)", "Birinapant|NSA (0.8)", "Birinapant|NSA (2)",
"Birinapant|NSA (5)", "Birinapant|NSA (12.5)", "Birinapant|NSA (0.8)|QVD-Oph (25)",
"Birinapant|NSA (2)|QVD-Oph (25)", "Birinapant|NSA (5)|QVD-Oph (25)", "Birinapant|NSA (12.5)|QVD-Oph (25)",
"Birinapant + Ipatasertib 2µM", "Birinpant + Ruxolitinib 2µM", "Birinapant + Bafilomycin A1 0.08µM",
"Birinapant + NSA 0.8µM",
"Birinapant + NSA 0.8µM + QVD-Oph 25µM", "Birinapant + NSA 2µM", "Birinapant + NSA 5µM",
"Birinapant + NSA 2µM + QVD-Oph 25µM", "Birinapant + NSA 5μM + QVD-Oph 25μM",
"DMSO", "empty", "Necrostatin-1|QVD-Oph")
Analysis
Overview - Stimulated
ld <- 1.5
cl3 <- 2.5
cl3small <- 2.25
apo <- 1.45
il2 <- 1.5
tnf <- 1.5
ifng <- 1.25
ki67 <- 1.5
gzmb <- 1.5
gmcsf <- 1.5
il10 <- 1.5
gzmb <- 1.5
plotTab <- ggcells(sce.stim, features = c("CD3", "FSCA", "SSCA", "Ki67", "Apotracker", "LD", "cleaved_caspase_3",
"IL2", "IL10", "TNF", "IFNg", "GMCSF"), exprs_values = "exprs")[[1]]
## Randomly sample to reduce overplotting
set.seed(100)
plotTab <- plotTab[sample(1:nrow(plotTab), nrow(plotTab)), ]
## Plot Treatment
p <- plotTab %>%
filter(!is.na(UMAP_FSCA.1)) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype %in% c("debris"), !is.na(celltype)) %>%
ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = Treatment)) +
geom_point(size = 0.1) +
xlab("UMAP1") + ylab("UMAP2") +
lgd + #scale_colour_viridis_c() +
theme_void() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none") +
scale_colour_manual(values = c("DMSO_stim" = "grey70", "Birinapant_low_stim" = "#AD1457"
))
p
ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_treatment_stim.png"),
height = 15, width = 15, units = "cm")
## Plot Diagnosis
p <- plotTab %>%
filter(!is.na(UMAP_FSCA.1)) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype %in% c("debris"), !is.na(celltype)) %>%
ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = Diagnosis_simple)) +
geom_point(size = 0.1) +
xlab("UMAP1") + ylab("UMAP2") +
lgd + #scale_colour_viridis_c() +
theme_void() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none") +
scale_colour_manual(values = c("T-LGL" = "#A5D6A7", "T-PLL" = "#64B5F6","healthy" = "#FFA000"
))
p
ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_diagnosis_stim.png"),
height = 15, width = 15, units = "cm")
## Plot celltype
p <- plotTab %>%
filter(!is.na(UMAP_FSCA.1)) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype %in% c("debris"), !is.na(celltype)) %>%
#mutate(celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = celltype_broad)) +
geom_point(size = 0.1) +
xlab("UMAP1") + ylab("UMAP2") +
lgd + #scale_colour_viridis_c() +
theme_void() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none") +
scale_colour_manual(values = c("T-LGL" = "#A5D6A7", "T-PLL" = "#64B5F6", "T-Cells" = "#0D47A1", "NK-Cells" = "#AB47BC",
"B-Cells" = "#FFA000"#, "monocyte" = "#66BB6A""#F44336"
))
p
ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_celltype_stim.png"),
height = 15, width = 15, units = "cm")
p1 <- plotTab %>%
filter(!is.na(UMAP_FSCA.1)) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype %in% c("debris"), !is.na(celltype)) %>%
mutate(celltype_broad = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
celltype_broad = ifelse(FSCA < 750000, "dead", celltype_broad),
celltype_broad = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3, "apoptotic", celltype_broad),
celltype_broad = ifelse(celltype == "dead" & cleaved_caspase_3 > cl3, "apoptotic", celltype_broad)) %>%
ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = celltype_broad)) +
geom_point(size = 0.1) +
xlab("UMAP1") + ylab("UMAP2") +
lgd +
theme_void() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none") +
scale_colour_manual(values = c("live" = "#A5D6A7","apoptotic" = "#FFA000",
"dead" = "#F44336"
))
p1
ggsave(plot = p1, filename = paste0(opt$plot, "TFlow_umap_celldeath_stim.png"),
height = 15, width = 15, units = "cm")
Overview - Unstimulated
plotTab <- ggcells(sce.unstim, features = c("CD3", "LD", "Apotracker", "cleaved_caspase_3", "FSCA", "CD25"), exprs_values = "exprs")[[1]]
## Randomly sample to reduce overplotting
set.seed(100)
plotTab <- plotTab[sample(1:nrow(plotTab), nrow(plotTab)), ]
ld <- 1.5
cl3 <- 2.5
cl3small <- 2.25
apo <- 1.45
## Plot cell types
p1 <- plotTab %>%
filter(!is.na(UMAP_FSCA.1)) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype %in% c("debris"), !is.na(celltype),
Treatment %in% c("DMSO", "Birinapant", "Selinexor", "QVDOph", "Birinapant|QVDOph")) %>%
ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = celltype_broad)) +
geom_point(size = 0.1) +
xlab("UMAP1") + ylab("UMAP2") +
lgd +
theme_void() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none") +
scale_colour_manual(values = c("T-LGL" = "#A5D6A7", "T-PLL" = "#64B5F6",
"T-Cells" = "#0D47A1", "NK-Cells" = "#AB47BC",
"B-Cells" = "#FFA000"
))
p1
ggsave(plot = p1, filename = paste0(opt$plot, "TFlow_umap_celltypes.png"),
height = 15, width = 15, units = "cm")
## Plot cell death
p1 <- plotTab %>%
filter(!is.na(UMAP_FSCA.1)) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype %in% c("debris"), !is.na(celltype),
Treatment %in% c("DMSO", "Birinapant", "Selinexor", "QVDOph", "Birinapant|QVDOph")) %>%
mutate(celltype_broad = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
celltype_broad = ifelse(FSCA < 750000, "dead", celltype_broad),
celltype_broad = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3, "apoptotic", celltype_broad),
celltype_broad = ifelse(celltype == "dead" & cleaved_caspase_3 > cl3, "apoptotic", celltype_broad)) %>%
ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = celltype_broad)) +
geom_point(size = 0.1) +
xlab("UMAP1") + ylab("UMAP2") +
lgd +
theme_void() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none") +
scale_colour_manual(values = c("live" = "#A5D6A7","apoptotic" = "#FFA000",
"dead" = "#F44336"
))
p1
ggsave(plot = p1, filename = paste0(opt$plot, "TFlow_umap_celldeath.png"),
height = 15, width = 15, units = "cm")
## Plot treatments
p <- plotTab %>%
filter(!is.na(UMAP_FSCA.1)) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype %in% c("debris"), !is.na(celltype),
Treatment %in% c("DMSO", "Birinapant", "Selinexor", "QVDOph", "Birinapant|QVDOph")) %>%
ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = Treatment)) +
geom_point(size = 0.05) +
xlab("UMAP1") + ylab("UMAP2") +
lgd +
theme_void() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none") +
scale_colour_manual(values = c("Birinapant" = "#AD1457", "QVDOph" = "#64B5F6",
"Birinapant|QVDOph" = "#43A047", "Selinexor" = "#FFD45F"
))
p
ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_treatments.png"),
height = 15, width = 15, units = "cm")
## Plot Diagnosis
p <- plotTab %>%
filter(!is.na(UMAP_FSCA.1)) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype %in% c("debris"), !is.na(celltype),
Treatment %in% c("DMSO", "Birinapant", "Selinexor", "QVDOph", "Birinapant|QVDOph")) %>%
ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = Diagnosis_simple)) +
geom_point(size = 0.1) +
xlab("UMAP1") + ylab("UMAP2") +
lgd + #scale_colour_viridis_c() +
theme_void() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none") +
scale_colour_manual(values = c("T-LGL" = "#A5D6A7", "T-PLL" = "#64B5F6","healthy" = "#FFA000"
))
p
ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_diagnosis.png"),
height = 15, width = 15, units = "cm")
Show cell death frequency relative to DMSO
plotTab <- ggcells(sce.unstim, features = c("CD3", "LD", "Apotracker", "cleaved_caspase_3", "FSCA", "CD25", "CD127"), exprs_values = "exprs")[[1]]
# Rare combinations of celltype x patient x treatment will be exclduded from the analysis, as this makes it more susceptible to noise.
countSub <- plotTab %>%
left_join(types, by = "uniqueBarcode") %>%
dplyr::count(patientID, Treatment, celltype_broad) %>%
filter(n < 100) %>%
mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype_broad)) %>%
pull(patTreatType) %>% unique()
plotTab %>%
select(patientID, Treatment) %>% unique() %>%
dplyr::count(Treatment)
## Treatment n
## 1 Birinapant 15
## 2 Birinapant|QVDOph 15
## 3 DMSO 15
## 4 QVDOph 15
## 5 Selinexor 15
compDrug <- c("Selinexor", "Birinapant_low", "Birinapant",
"Birinapant|QVDOph", "QVDOph",
"DMSO_stim", "Birinapant_low_stim")
ld <- 1.5
cl3 <- 2.5
cl3small <- 2.25
apo <- 1.45
testTab <- mclapply(unique(plotTab$patientID), mc.cores = ncores, function(x) {
testTab <- plotTab %>% filter(patientID == x) %>%
left_join(types, by = "uniqueBarcode") %>%
mutate(celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
mutate(cellDeath = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
cellDeath = ifelse(FSCA < 750000, "dead", cellDeath),
cellDeath = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3small, "apoptotic", cellDeath),
cellDeath = ifelse(cellDeath == "dead" & cleaved_caspase_3 > cl3, "apoptotic", cellDeath)) %>%
mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype_broad)) %>%
filter(!patTreatType %in% countSub) %>% # filter only for combinations with enough cells
group_by(patientID, celltype_broad, Treatment) %>%
mutate(nCells = length(uniqueBarcode)) %>%
ungroup() %>%
group_by(patientID, celltype_broad, cellDeath, Treatment) %>%
mutate(nCellsSub = length(uniqueBarcode), rel = nCellsSub / nCells) %>%
select(patientID, Diagnosis_simple, celltype_broad, Treatment, cellDeath, rel) %>% unique() %>%
pivot_wider(names_from = "cellDeath", values_from = "rel") %>%
mutate(Treatment = factor(Treatment, levels = c("DMSO", compDrug))) %>%
pivot_longer(cols = c("live", "apoptotic", "dead"),
names_to = "cellDeath", values_to = "rel") %>%
mutate(rel = ifelse(is.na(rel), 0, rel))
}) %>% bind_rows()
## Normalize to DMSO
viabTab <- mclapply(unique(testTab$patientID), mc.cores = ncores, function(x) {
lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
lapply(c("live", "apoptotic", "dead"), function(y) {
message(x)
message(z)
message(y)
testTab <- testTab %>% filter(patientID == x, cellDeath == y, celltype_broad == z)
dmso.val <- unique(testTab[testTab$Treatment == "DMSO", ]$rel)
if(length(dmso.val == 1)) {
testTab %>% mutate(viabNorm = rel / dmso.val,
viabDiff = rel - dmso.val)
} else {
testTab %>% mutate(viabNorm = NA)
}
}) %>% bind_rows()
}) %>% bind_rows()
}) %>% bind_rows()
viabTab.back <- viabTab
compDrug <- c("Selinexor", "Birinapant_low", "Birinapant",
"Birinapant|Q-VD-OPh", "Q-VD-OPh",
"DMSO_stim", "Birinapant_low_stim")
# Out of range observations
viabTab <- viabTab %>% mutate(viabNorm = ifelse(cellDeath == "live" & viabNorm > 1.25, 1.25, viabNorm),
viabNorm = ifelse(cellDeath == "apoptotic" & viabNorm > 5, 5, viabNorm),
viabNorm = ifelse(cellDeath == "dead" & viabNorm > 5, 5, viabNorm)) %>%
mutate(Treatment = str_replace(string = Treatment, pattern = "QVDOph", replacement = "Q-VD-OPh"))
viabTabOOR <- viabTab[viabTab$cellDeath == "live" & viabTab$viabNorm == 1.25 & viabTab$Treatment %in% c("Birinapant", "Selinexor", "Q-VD-OPh", "Birinapant|Q-VD-OPh"), ]
viabTabOOR <- viabTabOOR[!is.na(viabTabOOR$patientID), ]
## Plot normalized viability
p <- viabTab %>%
filter(cellDeath == "live", !is.na(viabNorm)) %>%
filter(Treatment %in% c("Selinexor", "Birinapant", "Q-VD-OPh", "Birinapant|Q-VD-OPh"
)) %>%
mutate(Treatment = factor(Treatment, levels = compDrug)) %>%
ggplot(aes(x = Treatment, y = viabNorm, fill = Treatment)) +
geom_boxplot(alpha = 0.6) +
geom_beeswarm(cex = 3.25, size = 2.5, shape = 21, data = viabTab[viabTab$viabNorm < 1.25 & viabTab$cellDeath == "live" & viabTab$Treatment %in% c("Selinexor", "Birinapant", "Q-VD-OPh", "Birinapant|Q-VD-OPh") & !is.na(viabTab$viabNorm), ]) +
geom_beeswarm(cex = 3, size = 2.25, shape = 24, data = viabTabOOR) +
#geom_text(aes(label = patientID)) +
geom_hline(yintercept = 1, linetype = "dashed") +
xlab("") + ylab("Living Cells") +
lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
legend.position = "none",
text = element_text(size = 20),
panel.border = element_rect(linewidth = 1, fill = "transparent")) +
facet_wrap(. ~ celltype_broad, nrow = 1) +
scale_fill_manual(values = c("Birinapant" = "#AD1457", "Selinexor" = "#FFD45F", "Q-VD-OPh" = "#64B5F6",
"Birinapant|Q-VD-OPh" = "#43A047"))
p
ggsave(plot = p, filename = paste0(opt$plot, "live_viabNorm_celltypes.png"),
height = 14, width = 30, units = "cm")
## Show the cell death type
viabTabOOR <- viabTab[viabTab$cellDeath == "apoptotic" & viabTab$viabNorm == 5 & viabTab$Treatment %in% c("Birinapant", "Selinexor", "Q-VD-OPh", "Birinapant|Q-VD-OPh"), ]
viabTabOOR <- viabTabOOR[!is.na(viabTabOOR$patientID), ]
## Apoptotic
p <- viabTab %>%
filter(cellDeath == "apoptotic",
!is.na(viabNorm)) %>%
filter(!Treatment %in% c("DMSO")) %>%
mutate(Treatment = factor(Treatment, levels = compDrug)) %>%
ggplot(aes(x = Treatment, y = viabNorm, fill = Treatment)) +
geom_boxplot(alpha = 0.6) +
geom_beeswarm(cex = 3.25, size = 2.5, shape = 21, data = viabTab[viabTab$viabNorm < 5 & viabTab$cellDeath == "apoptotic" & viabTab$Treatment %in% c("Selinexor", "Birinapant", "Q-VD-OPh", "Birinapant|Q-VD-OPh") & !is.na(viabTab$viabNorm), ]) +
geom_beeswarm(cex = 3, size = 2.25, shape = 24, data = viabTabOOR) +
geom_hline(yintercept = 1, linetype = "dashed") +
xlab("") + ylab("Apoptotic Cells") +
lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
legend.position = "none",
text = element_text(size = 20),
panel.border = element_rect(linewidth = 1, fill = "transparent")) +
facet_wrap(. ~ celltype_broad, nrow = 1) +
scale_fill_manual(values = c("Birinapant" = "#AD1457", "Selinexor" = "#FFD45F", "Q-VD-OPh" = "#64B5F6",
"Birinapant|Q-VD-OPh" = "#43A047"))
p
ggsave(plot = p, filename = paste0(opt$plot, "apoptotic_viabNorm_celltypes.png"),
height = 14, width = 30, units = "cm")
## Non-apoptotic
viabTabOOR <- viabTab[viabTab$cellDeath == "dead" & viabTab$viabNorm == 5 & viabTab$Treatment %in% c("Birinapant", "Selinexor", "Q-VD-OPh", "Birinapant|Q-VD-OPh"), ]
viabTabOOR <- viabTabOOR[!is.na(viabTabOOR$patientID), ]
p <- viabTab %>%
filter(cellDeath == "dead", #celltype_broad %in% c("lymphoma"),
!is.na(viabNorm)) %>%
filter(!Treatment %in% c("DMSO")) %>%
mutate(Treatment = factor(Treatment, levels = compDrug)) %>%
ggplot(aes(x = Treatment, y = viabNorm, fill = Treatment)) +
geom_boxplot(alpha = 0.6) +
geom_beeswarm(cex = 3.25, size = 2.5, shape = 21, data = viabTab[viabTab$viabNorm < 5 & viabTab$cellDeath == "dead" & viabTab$Treatment %in% c("Selinexor", "Birinapant", "Q-VD-OPh", "Birinapant|Q-VD-OPh") & !is.na(viabTab$viabNorm), ]) +
geom_beeswarm(cex = 3, size = 2.25, shape = 24, data = viabTabOOR) +
geom_hline(yintercept = 1, linetype = "dashed") +
xlab("") + ylab("Non-Apoptotic Dying Cells") +
lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
legend.position = "none",
text = element_text(size = 20),
panel.border = element_rect(linewidth = 1, fill = "transparent")) +
facet_wrap(. ~ celltype_broad, nrow = 1) +
scale_fill_manual(values = c("Birinapant" = "#AD1457", "Selinexor" = "#FFD45F", "Q-VD-OPh" = "#64B5F6",
"Birinapant|Q-VD-OPh" = "#43A047"))
p
ggsave(plot = p, filename = paste0(opt$plot, "non_apoptotic_viabNorm_celltypes.png"),
height = 14, width = 30, units = "cm")
Perform T-test
de.res <- lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
lapply(c("live", "apoptotic", "dead"), function(y) {
lapply(c("Birinapant", "Selinexor", "Q-VD-OPh", "Birinapant|Q-VD-OPh"
), function(x) {
message(z)
message(y)
message(x)
testTab <- viabTab %>% filter(celltype_broad == z, cellDeath == y, Treatment %in% c("DMSO", x)) %>%
mutate(Treatment = factor(Treatment, levels = c("DMSO", x))) %>%
select(-viabNorm, -viabDiff) %>%
pivot_wider(names_from = "Treatment", values_from = "rel")
if(nrow(testTab) >= 1) {
colnames(testTab)[[which(colnames(testTab) == x)]] <- "compTreat"
#testTab <- testTab[complete.cases(testTab), ]
res <- t.test(testTab$compTreat, testTab$DMSO, alternative = "two.sided",
var.equal = TRUE, paired = TRUE)
resTab <- data.frame(celltype_broad = z, cellDeath = y, p = res$p.value, diff = res$estimate, Treatment = x)
}
# t.test(testTab$compTreat, alternative = "less",
# var.equal = TRUE, mu = 1)
}) %>% bind_rows()
}) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"),
p.adj = round(p.adj, 3),
diff = round(diff, 2))
rownames(de.res) <- NULL
de.res
## celltype_broad cellDeath p diff Treatment p.adj
## 1 Lymphoma Cells live 2.127039e-03 -0.10 Birinapant 0.009
## 2 Lymphoma Cells live 1.180840e-02 -0.23 Selinexor 0.023
## 3 Lymphoma Cells live 1.074801e-01 0.06 Q-VD-OPh 0.129
## 4 Lymphoma Cells live 3.270611e-03 -0.44 Birinapant|Q-VD-OPh 0.010
## 5 Lymphoma Cells apoptotic 7.680785e-01 0.01 Birinapant 0.801
## 6 Lymphoma Cells apoptotic 2.414977e-02 0.12 Selinexor 0.040
## 7 Lymphoma Cells apoptotic 1.243085e-02 -0.12 Q-VD-OPh 0.023
## 8 Lymphoma Cells apoptotic 8.204419e-01 0.01 Birinapant|Q-VD-OPh 0.838
## 9 Lymphoma Cells dead 3.120977e-04 0.09 Birinapant 0.003
## 10 Lymphoma Cells dead 1.288734e-02 0.11 Selinexor 0.023
## 11 Lymphoma Cells dead 1.035210e-01 0.06 Q-VD-OPh 0.127
## 12 Lymphoma Cells dead 5.487479e-04 0.43 Birinapant|Q-VD-OPh 0.004
## 13 B-Cells live 1.138709e-02 -0.14 Birinapant 0.023
## 14 B-Cells live 1.194643e-03 -0.35 Selinexor 0.006
## 15 B-Cells live 2.115535e-02 0.04 Q-VD-OPh 0.036
## 16 B-Cells live 2.706574e-02 -0.20 Birinapant|Q-VD-OPh 0.042
## 17 B-Cells apoptotic 6.132582e-02 0.05 Birinapant 0.089
## 18 B-Cells apoptotic 4.441851e-03 0.30 Selinexor 0.013
## 19 B-Cells apoptotic 9.382079e-04 -0.17 Q-VD-OPh 0.006
## 20 B-Cells apoptotic 1.934377e-03 -0.13 Birinapant|Q-VD-OPh 0.008
## 21 B-Cells dead 4.671281e-02 0.09 Birinapant 0.070
## 22 B-Cells dead 9.250892e-02 0.05 Selinexor 0.123
## 23 B-Cells dead 9.290184e-03 0.13 Q-VD-OPh 0.019
## 24 B-Cells dead 6.522731e-03 0.33 Birinapant|Q-VD-OPh 0.016
## 25 T-Cells live 9.522270e-02 -0.08 Birinapant 0.124
## 26 T-Cells live 9.186048e-04 -0.22 Selinexor 0.006
## 27 T-Cells live 1.880645e-03 0.10 Q-VD-OPh 0.008
## 28 T-Cells live 3.326217e-01 -0.09 Birinapant|Q-VD-OPh 0.355
## 29 T-Cells apoptotic 1.007910e-01 0.04 Birinapant 0.127
## 30 T-Cells apoptotic 2.171900e-04 0.17 Selinexor 0.003
## 31 T-Cells apoptotic 7.737700e-05 -0.15 Q-VD-OPh 0.002
## 32 T-Cells apoptotic 7.070642e-03 -0.10 Birinapant|Q-VD-OPh 0.016
## 33 T-Cells dead 3.108288e-01 0.04 Birinapant 0.339
## 34 T-Cells dead 1.552538e-01 0.05 Selinexor 0.177
## 35 T-Cells dead 7.977198e-02 0.05 Q-VD-OPh 0.109
## 36 T-Cells dead 7.603422e-02 0.18 Birinapant|Q-VD-OPh 0.107
## 37 NK-Cells live 1.326006e-04 -0.08 Birinapant 0.002
## 38 NK-Cells live 5.257889e-03 -0.14 Selinexor 0.013
## 39 NK-Cells live 5.290138e-03 0.10 Q-VD-OPh 0.013
## 40 NK-Cells live 1.125779e-01 -0.07 Birinapant|Q-VD-OPh 0.132
## 41 NK-Cells apoptotic 3.057807e-03 0.06 Birinapant 0.010
## 42 NK-Cells apoptotic 1.221021e-05 0.10 Selinexor 0.001
## 43 NK-Cells apoptotic 2.714987e-03 -0.15 Q-VD-OPh 0.010
## 44 NK-Cells apoptotic 9.244809e-01 0.00 Birinapant|Q-VD-OPh 0.924
## 45 NK-Cells dead 2.703958e-02 0.02 Birinapant 0.042
## 46 NK-Cells dead 2.209660e-01 0.04 Selinexor 0.247
## 47 NK-Cells dead 5.101816e-03 0.06 Q-VD-OPh 0.013
## 48 NK-Cells dead 8.600919e-03 0.08 Birinapant|Q-VD-OPh 0.019
## Combination index
lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
lapply(c("live")[[1]], function(x) {
lapply(c("Birinapant|QVDOph"), function(y) {
message(z)
message(x)
message(y)
testTab <- viabTab.back %>%
filter(celltype_broad == z, cellDeath == x, Treatment %in% c("DMSO", "QVDOph", "Birinapant", y)) %>%
select(-rel, -viabDiff) %>%
pivot_wider(names_from = "Treatment", values_from = "viabNorm")
colnames(testTab)[[which(colnames(testTab) == y)]] <- "compTreat"
if(y == "Birinapant|QVDOph"){
testTab <- testTab %>% mutate(expViab = Birinapant * QVDOph)
}
## Remove samples with two NAs
testTab <- testTab %>% mutate(#compTreat = ifelse(is.na(compTreat), 0, compTreat),
CI = compTreat / expViab) %>% filter(!is.na(CI))
res <- t.test(testTab$expViab, testTab$CI, alternative = "two.sided", var.equal = TRUE, paired = TRUE)
resTab <- data.frame(celltype_broad = z, cellDeath = x, Treatment = y, diff = res$estimate, p = res$p.value,
meanCI = mean(testTab$CI))
}) %>% bind_rows()
}) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"),
p.adj = round(p.adj, 3),
meanCI = round(meanCI, 2))
## celltype_broad cellDeath Treatment diff
## mean difference...1 Lymphoma Cells live Birinapant|QVDOph 0.447124843
## mean difference...2 B-Cells live Birinapant|QVDOph -0.576507228
## mean difference...3 T-Cells live Birinapant|QVDOph 0.052968832
## mean difference...4 NK-Cells live Birinapant|QVDOph -0.008682567
## p meanCI p.adj
## mean difference...1 0.01631067 0.39 0.065
## mean difference...2 0.41684653 1.38 0.834
## mean difference...3 0.68971695 0.88 0.898
## mean difference...4 0.89829792 0.95 0.898
Show the proportions among live cells
countSub.sub <- plotTab %>%
left_join(types, by = "uniqueBarcode") %>%
mutate(celltype = ifelse(celltype %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype)) %>%
mutate( cellDeath = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
cellDeath = ifelse(FSCA < 750000, "dead", cellDeath),
cellDeath = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3small, "apoptotic", cellDeath),
cellDeath = ifelse(cellDeath == "dead" & cleaved_caspase_3 > cl3, "apoptotic", cellDeath)) %>%
filter(cellDeath == "live") %>%
dplyr::count(patientID, Treatment, celltype) %>%
filter(n < 100) %>%
mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype)) %>%
pull(patTreatType) %>% unique()
#### Now look at the proportions of healthy cells
testTab <- mclapply(unique(plotTab[plotTab$Diagnosis_simple == "healthy", ]$patientID), mc.cores = ncores, function(x) {
plotTab %>% filter(patientID == x) %>%
left_join(types, by = "uniqueBarcode") %>%
filter(!celltype_broad %in% c("debris")) %>%
mutate(celltype = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype),
celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
mutate(cellDeath = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
cellDeath = ifelse(FSCA < 750000, "dead", cellDeath),
cellDeath = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3small, "apoptotic", cellDeath),
cellDeath = ifelse(cellDeath == "dead" & cleaved_caspase_3 > cl3, "apoptotic", cellDeath)) %>%
filter(cellDeath == "live") %>%
mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype)) %>%
filter(!patTreatType %in% countSub.sub) %>% # filter only for combinations with enough cells
group_by(patientID, Treatment) %>%
mutate(nCells = length(uniqueBarcode)) %>%
ungroup() %>%
group_by(patientID, celltype, Treatment) %>%
mutate(nCellsSub = length(uniqueBarcode), rel = nCellsSub / nCells) %>%
select(patientID, Diagnosis_simple, celltype_broad, celltype, Treatment, nCellsSub, nCells, cellDeath, rel) %>% unique()
}) %>% bind_rows()
freqTab <- mclapply(unique(testTab$patientID), mc.cores = ncores, function(x) {
lapply(unique(testTab$celltype), function(z) {
lapply(c("live"), function(y) {
message(x)
message(z)
message(y)
testTab <- testTab %>% filter(patientID == x, cellDeath == y, celltype == z)
dmso.val <- unique(testTab[testTab$Treatment == "DMSO", ]$rel)
if(length(dmso.val == 1)) {
testTab %>% mutate(viabNorm = rel / dmso.val,
viabDiff = rel - dmso.val)
} else {
testTab %>% mutate(viabNorm = NA)
}
}) %>% bind_rows()
}) %>% bind_rows()
}) %>% bind_rows()
lapply(c("Birinapant", "Selinexor"), function(x) {
lapply(c("T-Cells"), function(y) {
df <- freqTab %>%
filter(cellDeath == "live", !celltype_broad %in% c("debris"),
Treatment == x, viabNorm != Inf) %>%
group_by(Treatment, celltype) %>%
summarise(meanViab = mean(viabNorm, na.rm = TRUE),
sdViab = sd(viabNorm, na.rm = TRUE) / sqrt(length(viabNorm))) %>%
ungroup()
cellOrder <- df %>%
arrange(desc(meanViab)) %>%
pull(celltype) %>%
unique()
p <- df %>%
mutate(celltype = factor(celltype, levels = cellOrder),
celltype_broad = ifelse(celltype %in% c("CD16- NK-cell", "CD16+ NK-cell"), "NK-Cells", NA),
celltype_broad = ifelse(celltype %in% c("naive B-cell", "memory B-cell"), "B-Cells", celltype_broad),
celltype_broad = ifelse(celltype %in% c("T-PLL", "T-LGL", "Lymphoma Cells"), "Lymphoma Cells", celltype_broad),
celltype_broad = ifelse(celltype %in% c("CD8 EM T-cell", "CD8 CM T-cell",
"CD8 T-emra", "naive CD8 T-cell", "gd T-cell", "CD4 CM T-cell",
"naive CD4 T-cell", "CD4 EM T-cell", "T-reg",
"CD4 T-emra"), "T-Cells", celltype_broad)) %>%
ggplot(aes(x = celltype, y = meanViab, fill = celltype_broad)) +
geom_bar(stat = "identity", colour = "black", alpha = 0.6) +
geom_beeswarm(aes(x = celltype, y = viabNorm), size = 3, shape = 21, cex = 2, data = freqTab[freqTab$Treatment == x, ]) +
geom_errorbar(aes(ymin = meanViab, ymax = meanViab + sdViab), width = 0.5) +
geom_hline(yintercept = 1, linetype = "dashed") +
xlab("") +
ylab("Frequency Rel. to DMSO") + ggtitle(paste0(x, " Effect on Healthy Donors (Living Cells)")) +
lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
legend.position = "none",
panel.border = element_rect(linewidth = 1, fill = "transparent"),
axis.text = element_text(size = 15),
text = element_text(size = 15)) +
scale_fill_manual(values = c("Lymphoma Cells" = "#64B5F6", "T-Cells" = "#0D47A1", "NK-Cells" = "#AB47BC",
"B-Cells" = "#FFA000"
))
ggsave(plot = p, filename = paste0(opt$plot, x, "_diff_effect_all_subtypes_bar.png"),
height = 14, width = 21, units = "cm")
p
})
})
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
## [[1]]
## [[1]][[1]]
##
##
## [[2]]
## [[2]][[1]]
Compute the fraction of proliferating cells
plotTab <- ggcells(sce.stim, features = c("FSCA", "Ki67", "Apotracker", "LD", "cleaved_caspase_3", "CD45RA", "CCR7",
"IL2", "IL10", "TNF", "IFNg", "GMCSF", "GranzymeB", "CD16", "CD56", "CD3", "CD19", "CD4", "CD8", "CD27", "TCRgd"), exprs_values = "exprs")[[1]]
## Get live cells
brc.live <- plotTab %>%
left_join(types, by = "uniqueBarcode") %>%
mutate(cellDeath = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
cellDeath = ifelse(FSCA < 750000, "dead", cellDeath),
cellDeath = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3small, "apoptotic", cellDeath),
cellDeath = ifelse(cellDeath == "dead" & cleaved_caspase_3 > cl3, "apoptotic", cellDeath)) %>%
filter(cellDeath %in% c("live")) %>%
pull(uniqueBarcode) %>% unique()
## Get only samples with decent cell counts
countSub <- plotTab %>%
filter(uniqueBarcode %in% brc.live) %>%
left_join(types, by = "uniqueBarcode") %>%
dplyr::count(patientID, Treatment, celltype_broad) %>%
filter(n < 100) %>%
mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype_broad)) %>%
pull(patTreatType) %>% unique()
compDrug <- c("Selinexor", "Birinapant_low", "Birinapant", "Birinapant|QVDOph", "QVDOph",
"DMSO_stim", "Birinapant_low_stim")
il2 <- 1.5
tnf <- 1.5
ifng <- 1.25
ki67 <- 1.5
gzmb <- 1.5
gmcsf <- 1.5
il10 <- 1.5
testTab <- mclapply(c("IL2", "Ki67", "TNF", "IFNg", "GMCSF", "IL10"), mc.cores = ncores,
function(x) {
if(x == "IL2") {cut.off <- il2}
if(x == "Ki67") {cut.off <- ki67}
if(x == "TNF") {cut.off <- tnf}
if(x == "IFNg") {cut.off <- ifng}
if(x == "GMCSF") {cut.off <- gmcsf}
if(x == "IL10") {cut.off <- il10}
plotTab %>%
left_join(types, by = "uniqueBarcode") %>%
mutate(celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
filter(FSCA < 3000000, uniqueBarcode %in% brc.live) %>%
pivot_longer(cols = x, names_to = "var", values_to = "expr") %>%
mutate(cellCytokine = ifelse(expr > cut.off, "pos", "neg")) %>%
mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype_broad)) %>%
filter(!patTreatType %in% countSub) %>% # filter only for combinations with enough cells
group_by(patientID, celltype_broad, Treatment) %>%
mutate(nCells = length(uniqueBarcode)) %>%
ungroup() %>%
group_by(patientID, celltype_broad, cellCytokine, Treatment) %>%
mutate(nCellsSub = length(uniqueBarcode), rel = nCellsSub / nCells) %>%
select(patientID, Diagnosis_simple, celltype_broad, Treatment, cellCytokine, rel) %>% unique() %>%
pivot_wider(names_from = "cellCytokine", values_from = "rel") %>%
mutate(Treatment = factor(Treatment, levels = c("DMSO", compDrug))) %>%
pivot_longer(cols = c("pos", "neg"),
names_to = "cellCytokine", values_to = "rel") %>%
mutate(rel = ifelse(is.na(rel), 0, rel)) %>%
mutate(cytokine = x)
}) %>% bind_rows()
## Normalize to DMSO
viabTab <- mclapply(c("IL2", "Ki67", "TNF", "IFNg", "GMCSF", "IL10"), mc.cores = ncores, function(y) {
lapply(unique(testTab$patientID), function(x) {
lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
message(x)
message(z)
#message(y)
testTab <- testTab %>% filter(patientID == x, cytokine == y, cellCytokine == "pos", celltype_broad == z)
dmso.val <- unique(testTab[testTab$Treatment == "DMSO_stim", ]$rel)
if(length(dmso.val == 1)) {
testTab %>% mutate(viabNorm = rel / dmso.val)
} else {
testTab %>% mutate(viabNorm = NA)
}
}) %>% bind_rows()
}) %>% bind_rows()
}) %>% bind_rows()
## Are there differences between stimulated and unstimulated DMSO
lapply(c("Ki67"), function(x) {
p <- viabTab %>% filter(cytokine == x) %>%
filter(Treatment %in% c("DMSO_stim", "Birinapant_low_stim")) %>%
mutate(Treatment = as.character(Treatment),
Treatment = ifelse(Treatment == "DMSO_stim", "DMSO", Treatment)) %>%
mutate(Treatment = ifelse(Treatment == "Birinapant_low_stim", "Birinapant", Treatment)) %>%
mutate(Treatment = factor(Treatment, levels = c("DMSO", compDrug))) %>%
ggplot(aes(x = Treatment, y = rel, fill = Treatment)) +
geom_line(linewidth = 0.5, aes(group = patientID)) +
geom_boxplot(alpha = 0.6) +
geom_beeswarm(cex = 3, size = 3, shape = 21) +
xlab("") + ylab(paste0("Frac. of ", x, " Expressing Cells")) +
ggtitle("Birinapant Effect on Proliferation") +
lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
legend.position = "none",
text = element_text(size = 17.5),
panel.border = element_rect(linewidth = 1, fill = "transparent"),
axis.text = element_text(size = 17.5)) +
scale_fill_manual(values = c("Birinapant" = "#AD1457", "DMSO" = "grey70")) +
facet_wrap(. ~ celltype_broad, nrow = 1)
if(x == "Ki67") {
ggsave(plot = p, filename = paste0(opt$plot, "TFlow_celltype_prol.png"),
height = 12, width = 30, units = "cm"
)
}
p
})
## [[1]]
Perform paired t-test
lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
lapply(c("IL2", "IL10", "TNF", "IFNg", "GMCSF", "Ki67")[[6]], function(y) {
lapply(c("Birinapant_low_stim"), function(x) {
message(z)
message(y)
message(x)
testTab <- viabTab %>% filter(celltype_broad == z, cellCytokine == "pos", cytokine == y, Treatment %in% c("DMSO_stim", x)) %>%
mutate(Treatment = factor(Treatment, levels = c("DMSO_stim", x))) %>%
select(-viabNorm) %>%
pivot_wider(names_from = "Treatment", values_from = "rel")
if(nrow(testTab) >= 1) {
colnames(testTab)[[which(colnames(testTab) == x)]] <- "compTreat"
#testTab <- testTab[complete.cases(testTab), ]
res <- t.test(testTab$compTreat, testTab$DMSO_stim, alternative = "two.sided",
var.equal = TRUE, paired = TRUE)
resTab <- data.frame(celltype_broad = z, cellDeath = y, p = res$p.value, diff = res$estimate, Treatment = x)
}
# t.test(testTab$compTreat, alternative = "less",
# var.equal = TRUE, mu = 1)
}) %>% bind_rows()
}) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"),
p.adj = round(p.adj, 3),
diff = round(diff, 2))
## celltype_broad cellDeath p diff
## mean difference...1 Lymphoma Cells Ki67 0.01894206 -0.02
## mean difference...2 B-Cells Ki67 0.91953751 0.00
## mean difference...3 T-Cells Ki67 0.09410205 -0.01
## mean difference...4 NK-Cells Ki67 0.52465156 0.00
## Treatment p.adj
## mean difference...1 Birinapant_low_stim 0.076
## mean difference...2 Birinapant_low_stim 0.920
## mean difference...3 Birinapant_low_stim 0.188
## mean difference...4 Birinapant_low_stim 0.700
lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
lapply(c("IL2", "TNF", "IFNg", "GMCSF"), function(y) {
lapply(c("Birinapant_low_stim"), function(x) {
message(z)
message(y)
message(x)
testTab <- viabTab %>% filter(celltype_broad == z, cellCytokine == "pos", cytokine == y, Treatment %in% c("DMSO_stim", x)) %>%
mutate(Treatment = factor(Treatment, levels = c("DMSO_stim", x))) %>%
select(-viabNorm) %>%
pivot_wider(names_from = "Treatment", values_from = "rel")
if(nrow(testTab) >= 1) {
colnames(testTab)[[which(colnames(testTab) == x)]] <- "compTreat"
#testTab <- testTab[complete.cases(testTab), ]
res <- t.test(testTab$compTreat, testTab$DMSO_stim, alternative = "two.sided",
var.equal = TRUE, paired = TRUE)
resTab <- data.frame(celltype_broad = z, cellDeath = y, p = res$p.value, diff = res$estimate, Treatment = x)
}
# t.test(testTab$compTreat, alternative = "less",
# var.equal = TRUE, mu = 1)
}) %>% bind_rows()
}) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"),
p.adj = round(p.adj, 3),
diff = round(diff, 2))
## celltype_broad cellDeath p diff
## mean difference...1 Lymphoma Cells IL2 0.340359218 0.01
## mean difference...2 Lymphoma Cells TNF 0.057969995 -0.03
## mean difference...3 Lymphoma Cells IFNg 0.350525204 -0.02
## mean difference...4 Lymphoma Cells GMCSF 0.022189361 0.01
## mean difference...5 B-Cells IL2 0.651648886 0.00
## mean difference...6 B-Cells TNF 0.004911755 -0.05
## mean difference...7 B-Cells IFNg 0.373888420 0.00
## mean difference...8 B-Cells GMCSF 0.912209574 0.00
## mean difference...9 T-Cells IL2 0.170823607 0.06
## mean difference...10 T-Cells TNF 0.252360329 -0.03
## mean difference...11 T-Cells IFNg 0.423812004 -0.01
## mean difference...12 T-Cells GMCSF 0.027517255 0.01
## mean difference...13 NK-Cells IL2 0.362919164 0.00
## mean difference...14 NK-Cells TNF 0.019890190 0.12
## mean difference...15 NK-Cells IFNg 0.103392260 0.13
## mean difference...16 NK-Cells GMCSF 0.016077270 0.17
## Treatment p.adj
## mean difference...1 Birinapant_low_stim 0.460
## mean difference...2 Birinapant_low_stim 0.155
## mean difference...3 Birinapant_low_stim 0.460
## mean difference...4 Birinapant_low_stim 0.088
## mean difference...5 Birinapant_low_stim 0.695
## mean difference...6 Birinapant_low_stim 0.079
## mean difference...7 Birinapant_low_stim 0.460
## mean difference...8 Birinapant_low_stim 0.912
## mean difference...9 Birinapant_low_stim 0.342
## mean difference...10 Birinapant_low_stim 0.449
## mean difference...11 Birinapant_low_stim 0.484
## mean difference...12 Birinapant_low_stim 0.088
## mean difference...13 Birinapant_low_stim 0.460
## mean difference...14 Birinapant_low_stim 0.088
## mean difference...15 Birinapant_low_stim 0.236
## mean difference...16 Birinapant_low_stim 0.088
Test for differential expression
sce.x <- sce.stim
compDrug <- c("Birinapant_low_stim")
colData(sce.x) <- colData(sce.x) %>%
data.frame() %>%
left_join(types, by = "uniqueBarcode") %>%
mutate(Treatment = factor(Treatment, levels = c("DMSO_stim", compDrug)),
cluster_id = celltype_broad,
celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
DataFrame()
testTreat <- c("Birinapant_low_stim")
cellSub <- c("Lymphoma Cells", "T-Cells", "B-Cells", "NK-Cells")
adtSub <- c("IL2", "GMCSF", "TNF", "IFNg")
## For some reason it still thinks there are NAs
sce.x <- sce.x[, !is.na(sce.x$celltype_broad)]
sce.x <- sce.x[, sce.x$uniqueBarcode %in% brc.live]
de.res <- lapply(testTreat, function(z) {
message(paste0("Fitting model for ", z))
lapply(cellSub, function(x) {
message(paste0("Celltype ", x))
sce.x <- sce.x[, sce.x$celltype_broad == x & sce.x$uniqueBarcode %in% brc.live]
## Set-up model formulas
meta.df <- unique(dplyr::select(data.frame(colData(sce.x)), patientID, Treatment, Stimulation))
design <- stats::model.matrix(~ patientID + Treatment, meta.df)
frml <- diffcyt::createFormula(meta.df, cols_fixed = c("patientID", "Treatment"))
## Get median counts
summed <- summarizeAssayByGroup(sce.x, ids = colData(sce.x)[,c("patientID", "Treatment")],
statistics = "median", assay.type = "exprs") # get raw counts and sum them up
colData(summed)$Treatment <- factor(colData(summed)$Treatment, levels = c("DMSO_stim", compDrug))
mat <- assay(summed, "median") %>% as.matrix()
colnames(mat) <- paste0(summed$patientID, ".", summed$Treatment)
## For every marker assemble the corresponding data frame and fit a linear model
lapply(adtSub, function(u) {
df <- data.frame(medianExpr = mat[u, ], patTreat = colnames(mat), ncells = summed$ncells) %>%
separate(col = "patTreat", into = c("patientID", "Treatment"), sep = "\\.")
testTab <- frml$data %>% left_join(df, by = c("patientID", "Treatment")) %>%
dplyr::rename(y = medianExpr) %>%
mutate(Treatment = factor(Treatment, levels = c("DMSO_stim", compDrug))) %>%
filter(!ncells < 100)
testTab
## Fit the linear model
fit <- lm(frml$formula, data = testTab, weights = testTab$ncells)
## Fit contrasts
cntr <- rep(0, times = length(names(fit$coefficients)))
cntr[[which(names(fit$coefficients) == paste0("Treatment", z))]] <- 1
cntr <- diffcyt::createContrast(cntr) %>% t()
res <- multcomp::glht(fit, cntr)
## Assemble output
resTab <- data.frame(celltype = x, Treatment = z, adt = u, estimate = summary(res)$test$coefficient, p = summary(res)$test$pvalues)
}) %>% bind_rows()
}) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"))
de.res %>%
mutate(estimate = round(estimate, 2),
p.adj = round(p.adj, 3))
## celltype Treatment adt estimate p p.adj
## 1...1 Lymphoma Cells Birinapant_low_stim IL2 0.10 0.278423532 0.318
## 1...2 Lymphoma Cells Birinapant_low_stim GMCSF 0.06 0.080219214 0.227
## 1...3 Lymphoma Cells Birinapant_low_stim TNF -0.10 0.130927389 0.246
## 1...4 Lymphoma Cells Birinapant_low_stim IFNg -0.04 0.318886545 0.340
## 1...5 T-Cells Birinapant_low_stim IL2 0.32 0.085254230 0.227
## 1...6 T-Cells Birinapant_low_stim GMCSF -0.03 0.258601059 0.318
## 1...7 T-Cells Birinapant_low_stim TNF -0.11 0.252182432 0.318
## 1...8 T-Cells Birinapant_low_stim IFNg -0.07 0.179300332 0.287
## 1...9 B-Cells Birinapant_low_stim IL2 0.00 0.983271082 0.983
## 1...10 B-Cells Birinapant_low_stim GMCSF 0.02 0.138553496 0.246
## 1...11 B-Cells Birinapant_low_stim TNF -0.12 0.004702416 0.073
## 1...12 B-Cells Birinapant_low_stim IFNg 0.02 0.039789897 0.159
## 1...13 NK-Cells Birinapant_low_stim IL2 -0.02 0.251544631 0.318
## 1...14 NK-Cells Birinapant_low_stim GMCSF 0.64 0.009137436 0.073
## 1...15 NK-Cells Birinapant_low_stim TNF 0.37 0.017227989 0.092
## 1...16 NK-Cells Birinapant_low_stim IFNg 0.23 0.129153616 0.246
p <- de.res %>%
mutate(pSign = -log10(p.adj) * sign(estimate),
label = ifelse(p.adj < 0.1, "*", "")) %>%
ggplot(aes(x = celltype, y = adt)) +
geom_tile(aes(fill = pSign), colour = "black") +
geom_text(aes(label = label), size = 6.5) +
xlab("") + ylab("") + ggtitle("Birinapant Effect on Cytokine\nExpression (< 10% FDR)") +
lgd + theme(panel.border = element_rect(colour = "black", fill=NA, linewidth = 1),
text = element_text(size = 17.5),
axis.text = element_text(size = 17.5),
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
scale_fill_gradient2(low = "#1976D2", high = "#F44336",
name = "-Log10(p.adj) \nwith Direction")
p
ggsave(plot = p, filename = paste0(opt$plot, "cytokine_heatmap.png"),
height = 15, width = 17.5, units = "cm")
Output session info
sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel grid stats4 stats graphics grDevices datasets
## [8] utils methods base
##
## other attached packages:
## [1] edgeR_4.0.16 limma_3.58.1
## [3] readxl_1.4.3 ggbeeswarm_0.7.2
## [5] circlize_0.4.16 ComplexHeatmap_2.18.0
## [7] gridExtra_2.3 BiocParallel_1.36.0
## [9] scater_1.30.1 scran_1.30.2
## [11] scuttle_1.12.0 SingleCellExperiment_1.24.0
## [13] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [15] GenomicRanges_1.54.1 GenomeInfoDb_1.38.8
## [17] IRanges_2.36.0 S4Vectors_0.40.2
## [19] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [21] matrixStats_1.3.0 lubridate_1.9.3
## [23] forcats_1.0.0 stringr_1.5.1
## [25] dplyr_1.1.4 purrr_1.0.2
## [27] readr_2.1.5 tidyr_1.3.1
## [29] tibble_3.2.1 ggplot2_3.5.1
## [31] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.2 bitops_1.0-7
## [3] cellranger_1.1.0 polyclip_1.10-7
## [5] XML_3.99-0.17 lifecycle_1.0.4
## [7] rstatix_0.7.2 doParallel_1.0.17
## [9] lattice_0.21-9 MASS_7.3-60
## [11] backports_1.5.0 magrittr_2.0.3
## [13] sass_0.4.9 rmarkdown_2.27
## [15] jquerylib_0.1.4 yaml_2.3.9
## [17] metapod_1.10.1 minqa_1.2.7
## [19] RColorBrewer_1.1-3 ConsensusClusterPlus_1.66.0
## [21] multcomp_1.4-25 abind_1.4-5
## [23] zlibbioc_1.48.2 Rtsne_0.17
## [25] RCurl_1.98-1.14 TH.data_1.1-2
## [27] tweenr_2.0.3 sandwich_3.1-0
## [29] GenomeInfoDbData_1.2.11 ggrepel_0.9.5
## [31] irlba_2.3.5.1 dqrng_0.4.1
## [33] DelayedMatrixStats_1.24.0 codetools_0.2-19
## [35] DelayedArray_0.28.0 ggforce_0.4.2
## [37] tidyselect_1.2.1 shape_1.4.6.1
## [39] farver_2.1.2 lme4_1.1-35.5
## [41] ScaledMatrix_1.10.0 viridis_0.6.5
## [43] jsonlite_1.8.8 GetoptLong_1.0.5
## [45] BiocNeighbors_1.20.2 survival_3.5-7
## [47] iterators_1.0.14 systemfonts_1.1.0
## [49] foreach_1.5.2 tools_4.3.2
## [51] ggnewscale_0.5.0 ragg_1.3.2
## [53] Rcpp_1.0.12 glue_1.7.0
## [55] SparseArray_1.2.4 xfun_0.45
## [57] withr_3.0.0 fastmap_1.2.0
## [59] boot_1.3-28.1 bluster_1.12.0
## [61] fansi_1.0.6 digest_0.6.36
## [63] rsvd_1.0.5 timechange_0.3.0
## [65] R6_2.5.1 textshaping_0.4.0
## [67] colorspace_2.1-0 utf8_1.2.4
## [69] generics_0.1.3 renv_1.0.7
## [71] S4Arrays_1.2.1 pkgconfig_2.0.3
## [73] gtable_0.3.5 RProtoBufLib_2.14.1
## [75] XVector_0.42.0 diffcyt_1.22.1
## [77] htmltools_0.5.8.1 carData_3.0-5
## [79] clue_0.3-65 scales_1.3.0
## [81] png_0.1-8 colorRamps_2.3.4
## [83] knitr_1.48 rstudioapi_0.16.0
## [85] reshape2_1.4.4 tzdb_0.4.0
## [87] rjson_0.2.21 nlme_3.1-163
## [89] nloptr_2.1.1 cachem_1.1.0
## [91] zoo_1.8-12 GlobalOptions_0.1.2
## [93] vipor_0.4.7 pillar_1.9.0
## [95] vctrs_0.6.5 ggpubr_0.6.0
## [97] BiocSingular_1.18.0 car_3.1-2
## [99] cytolib_2.14.1 beachmat_2.18.1
## [101] cluster_2.1.4 beeswarm_0.4.0
## [103] evaluate_0.24.0 mvtnorm_1.2-5
## [105] cli_3.6.3 locfit_1.5-9.10
## [107] compiler_4.3.2 rlang_1.1.4
## [109] crayon_1.5.3 ggsignif_0.6.4
## [111] labeling_0.4.3 FlowSOM_2.10.0
## [113] plyr_1.8.9 flowCore_2.14.2
## [115] stringi_1.8.4 viridisLite_0.4.2
## [117] munsell_0.5.1 Matrix_1.6-1.1
## [119] hms_1.1.3 sparseMatrixStats_1.14.0
## [121] statmod_1.5.0 highr_0.11
## [123] igraph_2.0.3 broom_1.0.6
## [125] bslib_0.7.0